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| ABSTRACT 

A method  for  the  aromatic  differentiation  of  computer  functions 
(subroutines)  written  in  a high  level  language  is  discussed. 

A theory  is  developed  to  show  that  most  functions  that  arise 
in  applications  can  be  differentiated  automatically.  It  is  shown  how 
one  can  take  a FORTRAN  function  (subroutine)  and,  with  the  aid  of  a 
precompiler,  obtain  a FORTRAN  subroutine  that  computes  the  original 
function  and  its  desired  derivatives. 

Implementation  of  two  types  of  differentiation  is  described  j iv  t 

1)  Automatic  Taylor  series  expansion  of  FORTRAN  programs.  ' 
< 2)  Automatic  Gradient  calculation  of  FORTRAN  functions. 
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AUTOMATIC  PI  FFERENTiATlOU  OI'  COMPUTER  PROGRAMS 

G,  Kedem 

Introdu  -non 

The  use  ot  re  • .rren  > red  if.o:.  . to  ompute  derivatives  is  not  a new 
idea.  We  can  trace  the.  idea  b 1 •.  j.  far  js  19  “>1  when  it  was  used  to  com- 
pute the  Emden  functions  by  means  of  Taylor  series  expansion  by  J.  R.  Airy 
(See  [1]).  This  idea  has  apparently  been  rediscovered  many  times.  In  1964 
R.  E.  Moore  [7]  showed  how  one  could  automatically  get  Taylor  series  expan- 
sions of  FORTRAN-like  expressions  to  solve  initial  value  problems.  See  [1,2, 
5,7,8,10).  The  automatic  computation  of  partial  derivatives  was  implemented 
in  1967  by  A.  Reiter  and  J.  Gray  [11,14]  and  later  by  J.  Wertz  [12],  D.  Kuba  and 
L.  B.  Rail  [13],  and  R.  E.  Pugh  [9].  These  are  programs  known  to  us  but  the 
list  is  probably  not  complete. 

This  paper  suggests  a way  to  extend  the  process  to  functions  that  can 

be  written  in  an  algebraic  computer  language  (FORTRAN,  ALGOL  and  so  on) 

namely:  piecewise  factorable  functions. 

The  theoretical  part  of  this  paper  was  written  because  we  heard  too  often 

the  statement:  "Oh,  I believe  you  can  differentiate  arbitrary  expressions,  but 

FORTRAN  programs  ?!".  We  then  describe  the  implementation  of  two  types  of 

differentiation:  Taylor  series  expansion  (TAYLOR)  and  gradient  computation 

(GRADIENT),  via  the  use  of  the  AUGMENT  precompiler  [3]. 

The  method  described  has  a few  distinct  advantages: 

1)  AUGMENT  is  a highly  portab'e  precompiler,  therefore  with  little 

work  it  can  be  implemented  on  almost  any  system. 

Sponsored  by  the  United  States  Army  under  Contract  No.  DAAG29 -7 5-C -0024 . 


2)  The  packages  give  a very  easy  way  to  interface  with  any  FORTRAN 


program. 

3)  Any  other  type  of  differentiation  can  be  implemented  easily  (about 

two  weeks  of  work). 

4)  Functions  and  subroutines  can  be  "differentiated''  separately, 

incorporated  into  larger  systems  and  used  as  part  of  the  library. 

1.1.  Theory. 

Let  us  look  at  computer  programs  for  the  evaluation  of  numerical 
functions  that  arise  in  applications. 

We  use  the  FORTRAN  language  but  the  following  discussion 
applies  to  any  other  high  level  algebraic  language.  We  ignore  the  fact 
that  computers  don’t  really  work  with  real  numbers  and  that  library  functions 
like  SIN  and  COS  compute  only  approximations  to  functions  we  have  in  mind. 

We  look  at  FORTRAN  subroutines  and  functions  that  evaluate  some 
mathematical  functions.  We  assume  that  no  I/O  is  involved  and  that  "random 
numbers  are  not  used.  All  such  routines  have  a few  features  in  common: 

1)  For  every  set  of  values  of  the  formal  arguments  there  is  a fixed, 
finite  sequence  of  instructions  executed  (provided  that  these  values  are  with- 
in the  domain  of  definition  of  the  function  and  we  use  a correct  subroutine). 

2)  If  we  regard  DO  loops,  logical  statements  and  GOTO  statements  only 
as  a convenient  tool  for  defining  sequences  of  instructions,  we  see  that  all 
such  sequences  consist  of  arithmetic  operations  and  calls  to  library  functions. 

3)  At  each  step  we  use  only  previously  defined  values. 

4)  Most  of  the  functions  we  compute  are  piecewise  differentiable  or 
actually  piecewise  analytic. 
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In  order  to  make  the  analysis  more  precise  we  use  an  abstract  model 
for  algebraic  computer  languages,  namely  the  factorable  functions.  But  we 
will  try  to  point  out  the  analogy  between  the  model  and  computer  programs 
and  what  the  statements  about  factorable  functions  mean  when  we  translate 
them  to  facts  about  computer  programs. 

We  will  use  subscripts  to  denote  sequences  and  superscripts  to  de- 
note components  of  a vector,  f.  is  the  ith  element  in  a sequence  and  tJ 
is  the  jth  component  of  a vector  t . 

1.2.  Definitions. 

Let  f be  a finite  set  of  real  functions  of  one  or  more  real  arguments 

including  the  identity  function.  We  call  i the  set  of  basic  library  functions 

Let  f be  a map  from  Rn  to  Rm.  We  call  f factorable  function 

if  and  only  if  there  exists  a finite  sequence  of  functions 

f f : D C Rn—  R that  satisfy  the  following  conditions: 

1 k — 

1)  fj(x)  = x1 , f 2(x)  = x2,  . . . , Mx)  = xn 

2)  fk-m+l(X)  = fl(xl>  •••>  ^ fk(x)  = fnl(x1’  xD) 

3)  for  every  i,n<i<k,  either  3 ge  I 3 

f.(x)  = g(f  (x),  . . . , f (x) ) < i, 

i jj  3s  L 

or 

f (x)  = C € R , 

we  call  such  a sequence  a basic  sequence  and  we  say  that  the  sequence 

f ,f,  is  a basic  representation  for  f . 

1 k 
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Remarks . 


])  We  will  not  use  the  properties  of  the  real  numbers  in  most  of  the  discus- 
sion to  follow,  so  one  could  apply  the  theory  to  functions  over  any  set. 

2)  The  analogy  to  FORTRAN  is  clear.  The  sequence  (f. } that  is  a basic 
representation  for  f corresponds  to  subroutine  that  computes  f.  The  first 
n places  correspond  to  the  arguments  and  the  last  m to  the  result.  The 
middle  part  corresponds  to  the  body  of  the  subroutine.  A basic  representation 
for  f can  actually  be  written  down  by  following  the  path  of  execution 
(assuming  no  IF  and  GOTO  statements  are  involved). 

3)  Many  different  sequences  can  represent  the  same  function. 

4)  It  is  easy  (but  important)  to  verify  that  for  every  i.n  < i < k,  the  sequence 
f j , . . . , f.  represents  f.  . 

5)  It  is  always  assumed  that  the  compositions  are  well  defined,  that  is: 

If  1 - g(fj  , . . . , f . j then  for  every  x in  the  domain 
is  — 

of  f ( f (x)  , f (x)  . . . . f.  (xj)  is  in  the  domain  of  g.  Only  if  this 
li~  h~  ’s 

condition  is  satisfied,  can  we  call  f factorable. 

6)  Functions  which  are  identically  constant  are  thought  of  as  functions  of 


many  arguments.  The  number  of  arguments  depends  on  the  context. 


Example: 

Here  is  an  example  of  a factorable  function.  Let  c be  the  set 
{+,  /t  sin,  cos,  exp}  and  f(x)  = cos(x)  * exp(sin(x))  + sin(x)**2. 

A sequence  which  is  a basic  representation  for  f is: 
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1 ) f (x)  X 


2)  l^[x)  - cosif  | (x i ) 

3)  f ^ C x ) - sin(f  j (x) ) 


4)  f,(x)  = exp(f  (x)) 

4 3 

5)  f _ (x)  = f,(x)  f (xl 

d 3 3 

6)  f,(x)  = f (xl  f (x) 

6 2 4 

7)  f (x)  = f,  (x)  + f (x) 

/ o b 

1.3.  Composition. 

n k K m 

Let  q:R‘  -♦  R and  h:R  -*•  R be  factorable  functions.  Let 
’qL1+k  and  hl’h2-'-hL2 

be  basic  representations  for  g and  h.  From  these  two  sequences  we  con- 
struct a third  sequence  that  we  call  the  composition  of  q ,...,q  with 

1 Lj4k 


hj, . . . , hj.  The  sequence  is  of  the  form  f^,  f^,  . . . , f 
for  i = 1, . . . , L + k,  and 


Ll+L2 


where  f.  = q . 


c if  h = c 

l 


c c R 


g(f  ,f  ,...,f  ) if  h = g(h.  ,...,h  ) 

Ll+jl  Vj2  VV  1 V V 


for  i = k+1,  k+2, . . . , L 


Remark: 


The  composition  of  q,,  . . . , q with  h , . . . , h is  simply 

1 Ll+k  1 L? 

overlaying  the  first  k terms  in  the  sequence  h . . . , h on  top  of  the 

1 2 

last  k terms  of  the  sequence  q ^ , . . . , qL 
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Lemma  1:  Let  q and  h be  factorable  functions  q:  Rn  — Rk;  h:  Rk  -*  Rm 

•et  qr  ■ ’ ' ,cil  +k  and  hj,  . . . , be  basic  representations  for  q and  h,  then 
f - h( q>  is  a factorable  function  and  the  composition  of  q 

1 ’ ’ ’ ’ ’ L+k 

b ,»•••,  h is  a basic  representation  for  f . 

2 

Proof:  The  proof  follows  indirectly  from  the  definitions. 


with 


Corollary  1 : The  set  of  factorable  functions  is  closed  under  finite  number  of 

substitutions. 

Corollary  2:  Let  3 be  a collection  of  factorable  functions.  Let  f^.  . . ,f 

be  a sequence  of  functions  of  the  following  form. 

1) 

2) 


r 1 , n 

fl  = x ’ ' • ■ > f„  = * 


for  n < i < L f.  is  of  one  of  the  forms 

a)  f.  = const. 

l 

b)  = g(fi  , . . . , f ) for  some  qclf  <i, 

J1  ‘s  is 

c)  f = q(f  , . . . ,f  ) for  some  q€  3 j , . . . , j < i . 

J1  Jt  1 l 

Then,  for  n < j < L,  f^  is  a factorable  function.  We  say  that  the  sequence 

f, , . . . ,f.  is  a factorable  representation  for  f . 

J > ) 


Remark. 


The  substitution  of  one  sequence  into  another  corresponds  to  a call 

to  subroutine  or  function  in  FORTRAN  and  the  sequence  f , . . . ,f  in 

1 L 

Corollary  2 corresponds  to  a subroutine  that  uses  other  subroutines  in  the 
computations. 
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i . 4 . Differentiation: 


r 


In  this  section  we  show  how  one  can  compute  total  or  partial  derivatives 
of  factorable  functions.  We  show  how  one  can  get  from  a basic  representation 
for  a factorable  function  to  a new  sequence  of  functions  which  is  a factorable 
representation  of  the  original  function  and  its  total  (or  partial  derivatives  by 
means  of  simple  replacement  operations. 

We  start  with  an  example: 

Let  Z and  f(x)  be  as  in  section  1.2:  that  is,  £=  {+,  * , /,  sin, 

cos,  exp}  and 

f(x)  = cos(x)  * exp(sin(x))  + sin(x)**2 


before 

the  following  : 

1) 

= I (I  i 

2) 

f2  = cos(  fj) 

3) 

f3  = sin(fj) 

4) 

f^  = exp(f3) 

5) 

f5  = f3  * f3 

6) 

f6  = f2  * f4 

7) 

f7  = f6  + f 5 
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We  replace  this  sequence  by 


l ) P,  (X, 

l 


F[(X,X’)=X'  ) 


Y(X,X")  £ R 


2) 

f2 

•osiFj)  , 

f2 

= -sin(F^l 

» * f;  > 

3) 

''  F3 

sinlF^I  , 

F3 

- 30S( 

Fl> 

* F|  ) 

4) 

' F4 

exp(  F ^ ) 

F4 

= F4* 

F' 

‘j 

) 

5 ) 

< Fr 

--  F * F,  . , 

F^ 

II 

TJ 

-M. 

F’ 

+ F1  * F 

5 

3 3 

5 

3 

3 

3 3 

6 ) 

< F, 

- F * F 

F' 

= F , * 

F' 

+ F'  * F 

6 

2 4 

6 

2 

4 

2 4 

7) 

' F, 

7 

= F.  + F , 

D 3 

F7 

= F6  + 

F' 

s 

> 

Please  Note: 


2 „2 


a)  The  new  sequence  is  a factorable  representation  for  a function  F ;R  — r 

bi  There  is  a "simple"  correspondence  between  the  original  sequence  and 

the  new  sequence.  Each  term  in  the  original  sequence  of  the  form  g(fj) 

2 2 

was  replaced  by  a term  G(F^)  where:  G:  R — R is  a factorable  function, 

and  F is  the  jth  term  in  the  new  sequence.  Similarly  terms 

j 

of  the  form  g(f.  , f.  ) were  replaced  by  terms  of  the  form  G(F.  ,F.  ) 


J1  ]2 


where  G : R 


fit  ) 

c)  Finally,  F(t  ,1)  = / 

T~f(t)  I 
dt  't=t 


J1  J2 

a factorable  function  which  corresponds  to  g . 

"V 

. or  in  general:  If  h is  a function  of  t, 


h(t0)  - h„  , - hit)  |t=t<)  = h'0  then 


The  following  discussion  describes  the  general  case. 


Let  f be  a set  of  basic  library  functions.  Let  T be  an  operator 
that  maps  functions  to  functions.  We  assume  the  following: 

1)  T is  defined  for  all  factorable  functions  and  if  f : D C Rn  — R then 
T[  f]  : E C Rn  ^ — R^'  ^ (k  is  the  degree  of  T).  E - D X R m. 


3)  t[  I 1 = I (I  is  the  identity  map  in  R3). 

1 k j 

4)  For  every  gt  £ (g  : D^R  -*■  R)  3 a factorable  function 

s • k k / k - 1 ) * q 

G:  EC  R — R , (E  = D X R ),  3 for  every  factorable  function 

f : DCRn-RS,  f (D)  C_  D 

T[  g(f 1 , fS)]  - Gaff1],  ...,  T[fS])  . 

[ 
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1 


5)  T[  const.]  is  a constant  vector  function  (T[  const.]  c R ) and  £ a 


factorable  function  C : R — R y Vet  R T[  c]  = C(c). 


Theorem  1 . 


Let  f and  T satisfy  the  above  conditions.  Let  f be  a factorable 
function  f :R " - R,  and  let  f , . . . , f^  be  a basic  representation  for  f,  then: 
T[f]  is  a factorable  function. 


moreover  if  we  replace  the  terms  f, , . . . , fT  by  F. , . . . , F where 


i)  for  1 < i < n 
ii)  for  n < i < L 


r • • • *‘i  * r • • • 7‘l 

F^  is  the  sequence  y}>^} 


i,  k 


r 


G,'V'-’V  l£  f‘s9‘<fv"V 


s 


C(c) 


if 


fi  = c 


then  we  get  a factorable  representation  for  T[  f]  . 


Proof: 


By 


induction  on  the  length  of  the  sequence  representing  f . If 


f:Rn  - R 


and  f , . . . , fT  represents  f and  L - n+1  then  either 
1 ' L 


i J s 

a)  f = g(x  , . . . , x ) for  some  g t £ 


or 


b)  f = const. 


incase  a we  replace  the  sequence  x^xi 2,...^  , g(x  ,...,x  ) by 


1,1  l,k 

x ’ , . . . ,x 


!1 , 1 lljrv  /-i/  .. 

,. . • , x ’ ,.  . . ,X  , G(x 


n,  k 


, . . . ,X  , . . . ,X  ) 


^The  function  G.  is  the  factorable  function  corresponding  to  the  basic 


library  function  g.. 
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where  G is  the  fa  torable  function  corresponding  to  g . By  Assumptions  3 


and  4 


3.  J 

Tjf | T|g(x  ,.  . . ,x  )]  = G(x 


jl’* 


Jl'k 


,1 


md  therefore  we  have  a factorable  representation  for  T[  f]  . In  case  b if  f r 

then  we  replace  , . . . , x , c by  x ’ , . . . , x ' , C(  c)  and  again  by 

Assumptions  Tf  c]  = C(c).  Let  f:Rn-R,  f j fL  represent  f,  L=n+J, 

assume  the  theorem  is  true  for  all  sequences  of  length  L = n+i,  i<  j.  f=f^<x  > • 

is  either  of  the  form  gT(f.  ,...,f.  ),  gT  « £ l,,.*-,!  <L  or  “const. 

L jl  ]s  L 

For  f = const  the  same  argument  as  for  L = n+1  applies.  Let  ij,...,F^ 

be  the  sequence  that  corresponds  to  f^ , . . . , f^.  f-f^  = Q(  f ^ , . . . , f l , so 

by  construction.  F = GIF 
L 


By  the  induction  hypothesis,  i 


Iff 


3. 


s J1 

T If  1 and  F 1 < i < s are  factorable,  therefore!  By 
3 1 ' “ “ 


3; 


Lemma  1,  GIF  , ...,F  ) is  factorable . By  Assumption  4, 

J1  Js 


T[fJ  - T[  g(  f ,...,f  )]  = G(T[f  J., ...  , T[f . ]) 
ji  Js  J1  Js 


= GIF.  ,.. . ,F  ) 
'1  Js 


n f.  n. 


As  an  immediate  consequence  of  this  theorem  we  have: 

Corollary  3:  The  theorem  is  true  for  factorable  functions  from  Rn  to  Rm  . 
We  also  have  the  following  extension. 

Theorem  2 : Let  T and  £ satisfy  the  above.  Let  7 be  a set  of  factorable 

functions.  Let  f be  a factorable  function  f:Rn  — R,  and  f^,.  . . ,f  a 
factorable  representation  for  f,  that  is:  f^  = x t . . . , fR  = x and  for 

n < i < L f,  is  of  one  of  the  forms 

— l 


and 

n 

. . ,x  ) 
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C)  t const. 

If  we  replace  this  sequence  by  the  sequence  F^,.  . . , where, 
F x ^ and  for  n < i <_  L 


G(F  , . . . ,F  ) if  f g(f  ) , g<  £ 

J1  Js  1 J1  Js 

F.  < T[q ]( F . ,F.  ) if  f = q(f  q<  3 

1 J1  1 J1  Jt 

C(c)  if  f s c 

then,  for  n<  f < L,  the  sequence  Fj,...,Ff  is  a factorable  representation  for 
T[  f ] . The  proof  of  this  theorem  is  the  same  as  the  previous  one.  We  only 
have  to  use  the  following  lemma. 

Lemma  3: 

Let  f be  a factorable  function  f : DC  R -►  R.  Then,  for  every 

.s  /s  rn  p a a 

factorable  function  f : DC  R -•■  R with  f(D)  _ D , 

T[  f (f1 , . . . , f")]  = T[  f]  (T[  f1]  , . . . , T[  fn]  ) . 

Proof: 

By  induction  on  L,  the  length  of  the  sequence  representing  f . 

If  L = n + 1 then  either  f = c or 
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f - g(x 


,x  ) 


i i ? . 


If 

is 

if 


f : c there  is  nothing  to  prove  and  if  f 
true  by  Assumptions  3 and  4 

Assume  the  lemma  is  true  for  all  L 
f ~ c,  again  there  is  nothing  to  prove  an 


g(x 


then  the  lernma 


n+j  j < i and  L.  = n+i,  then 

d if  f = g(f  , ...  ,f  ) 

J 1 J o * 


j < L g c S'  then: 
s 

^1  ^ 1 /Nn  ^1 

T[f(f  , . . . , f )]  = T[g(  f (f  ),...,f  (f  ,...,f  ))] 

J1  s 

1 /\  /s  j ^ n 

= G(T[f  (f  ,.••  ,f  )],...  ,T[f.  (f  ,---»f  >1)  by  Assumption  4 
J1  Js 

- G(T[f  ] ( T [ f 1 ], . . . , T [ f ' 1 ] ) , . . . ,T[f.  l(T[f  1 1 , . . . ,T[fn]))  by  induction 

J1  •’s  hypotheses 

= T[g(f  , ...,f  )](T[f [],  . . . ,T[f  nD  by  Assumption  4 

J1  Js 

= T[f](T[f X], . . . , T [ f n ] ) =>  Q.  E.  D. 


Remarks : 

1)  Take  r to  be  the  set  of  standard  Fortran  functions  and  the  arithmetic 
operations.  Let  T be  an  operator  we  wish  to  implement  and  assume  that 
Assumptions  1-5  are  satisfied.  Theorem  1 then  says  the  following:  Let  F 
be  a subroutine  that  computes  a function  f (For  the  moment  we  assume  that 
no  IF  or  GOTO  statements  are  used  ).  Suppose  we  replace  each  variable 
by  a K vector  (L  vector  by  a KX  L-array  and  so  on),  replace  each  constant 
by  a constant  vector,  and  replace  each  arithmetic  operation  and  function  call 
by  a corresponding  subroutine.  Then  the  result  would  be  a Fortran  subroutine  that 
computes  T [ f ] . The  replacement  rules  are  simple  and  can  be  carried  out 
mechanically  by  a precompiler. 


j 
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<d)  Theorem  Z says  that  if  we  are  using  other  subroutines  or  functions 
in  the  course  of  the  computation  we  can  transform  them  separately.  In  the 
replacement  process  we  replace  the  original  call  by  a call  to  the  corresponding 
transformed  routine. 

1.5  Example  1.  Taylor  Series  Expansion 

1 d1 

Let  (X(t  )].  denote  — — — X(t)|  ' Assume  we  know  the  vdlues 

0 1 1 ' dt1  1 0 

of  lX(tQ)].,  [Y(tQ) ).,  J = 0,...,k. 

i)  If  Z = X ± Y , then  we  can  compute  [Z(tQ)].  by 

iHtylj  = [X(t0)].  ± [Y(t0)]j  , j =0,1,. ..,k» 

ii)  If  Z - X * Y , then  by  Leibnitz'  rule 

= t [X(,0»)s  ■ lY|t0»li-s  ’ 1 s0"--’ k ■ 

s=0 

iii)  If  Z = X/Y  and  Y(t  )#0,  then 

1-1 

[Z]  = l/Y(tQ)  • { [X(  tQ ) ] j XtZ<t0)]s-  [Y(t0)]._s}  j = 1, 2, . . . , k . 

J s =0 

iv)  If  Z = EXP(X),  then  we  can  compute  lZ(tQ)].  using  the  recursion  relations 

j-1 

[Z(t0)J  = Z «i-s)/j)[Z(t0)]s  • [X(t0>]  j = 1, 2, . . . , k . 

s=0 

It  is  well  known  that  there  are  recurrsion  relations  that  enable  one  to  com- 
pute successive  derivatives  of  functions  that  satisfy  rational  differential 
equations.  In  Appendix  A we  give  a list  of  such  recurrsion  relations  for  the 
most  common  special  functions.  (See  also  (1,7]). 

As  a matter  of  fact,  the  discussions  in  [1]  and  [7]  show  that  one  can 
write  such  recurrsion  relations  for  functions  satisfying  differential  equations 
of  the  form  Y’  = f( t,  Y)  where  f is  a factorable  function  and  £ is  a set 
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containing  the  arithmetic  operations  and  functions  satisfying  rational  differ- 


ential equations. 

Let  S*  be  a set  that  consists  of  the  identity  function,  the  arithmetic 
operations,  functions  that  satisfy  first  order  rational  differential  equations 
(like  exp  and  log),  pairs  of  functions  satisfying  second  order  rational  differ- 
ential equations  (like  sin  and  cos),  and  so  on.  Choose  the  set  <*  in  such  a 
way  that  all  recurrsion  relations  between  successive  derivatives  can  be 
expressed  as  factorable  functions.  For  example:  the  set  of  functions  in 
Appendix  A is  such  a set. 

s 

Let  k > 1 be  an  integer.  For  each  basic  library  function  g:  D C R — R 
take  G to  be  the  factorable  function  satisfying: 
a)  G : E = D X RS'  k -*  Rk+1 


s k 

b)  for  every  function  X : R -*  R , XeC  , if  X(  ) t D and 


[X(t0)].  = X.  « R , i = 0,1,...  ,k 


G(XQ,Xl,..  . ,Xk) 


i9<x<y»0 

[9<X(t  »), 


ls<x<t0»Ik 


We  define  T,  as  follows: 
k 


s k 

Let  f:  W C R -*  R be  of  class  C . 

Let  T[f]  = F be  the  function  satisfying: 

. r-»  ^ r,S-  k „k  + l 

a)  F : W X R -*  R 

s k 

b)  For  every  function  X : R -►  R : X t C , and  for  every  t R 
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i 


such  that  X< tQ)  t W 

Z [«xv'o  ^ 

| l(IXIto'»l,  \ 

n|x(«0,]o,  IX(.0»], I»VV  I : 

V [f<xV»k 

It  is  not  hard  to  see  that  the  operator  T and  the  set  P defined  aDove 
satisfy  assumptions  1-5  of  Theorem  1. 

Note  that  the  set  Z contains  only  analytic  functions  and  therefore 
the  factorable  functions  are  analytic.  In  the  next  section  we  will  show  that 
it  is  possible  to  include  in  the  basic  library  set,  functions  that  are  only 
piecewise  analytic  (like  max,  min  and  abs). 

Also  note  that  in  order  to  compute  high  order  derivatives  of  a func- 
tion one  has  to  use  all  the  lower  order  derivatives  of  that  function.  There- 

(k) 

fore  the  operator  T : f — f , k>  1 does  not  satisfy  assumptions  1-5  of 
Theorem  1. 

1.6  Example  2:  Partial  derivatives 

9f 

Let  f,.,  denote  

(J)  ax. 

Let  p be  as  in  Example  1. 

Let  k > 1 be  an  integer. 

s 

For  each  basic  library  function  g:  D (_  R -*  R take  G to  be  the 
factorable  function  satisfying: 

a)  G : D X RS‘  k - Rk+1 
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k s 

b)  For  every  function  h:  R — R , hi  C if  h(X^ ) i D then 


/glhix0l)  ' 


G(h(xo>  • W *Vxo»  - 


We  define  T,  as  follows: 
k 

s 1 

Let  f:  W C R -*  R be  of  class  C , let  T^[f]  = F be  the  function 
satisfying  : 

a)  F : W X RS  ‘ k - Rk+1 

k s 1 k 

b)  For  every  function  h:  R — R h e C , and  for  every  « R 

such  that  h(XQ ) t W 


io"i'(iro"-'-,‘,(kro 


Again  it  is  not  hard  to  see  that  the  above  4"  and  satisfy  the  assump- 


tions of  Theorem  1. 


m 

1.7.  Piecewise  factorable  functions: 

So  far  the  model  does  not  describe  computer  programs  which 
include  IF  and  GOTO  statements.  The  extension  is  straightforward. 

We  use  the  following  definition: 

Definition: 

A function  f : D C Rn  - Rm  is  a piecewise  factorable  function 

k 

if  and  only  if  there  exists  finite  number  of  sets  U. , . . . , U 3 D C U U 

J=1  1 

and  f restricted  to  each  is  a factorable  function.  A very  useful  fact 
about  piecewise  factorable  functions  is: 

Lemma  4: 

Let  q and  h be  piecewise  factorable  functions,  q : CRn»  Rk, 

h : D2C  Rk  - Rm  if  q(Dj)  = (ze  Rk  | z = q(x),  Xf  Dj}  C D2  then 
f : Dj  C Rn  - Rm  defined  by  f(x)  = h(q(x)),  x«  is  a piecewise  factor- 
able function. 

Proof: 

s 

Let  U . . . , U C Rn  be  sets  9 D C U U and  q I ..  is  a 
1 s~  k 1 i=i  1 t '• 

factorable  function.  Let  tL ..... V CR'  be  sets  $ L*.  V. 

1 1 j=  1 1 

and  h | is  a factorable  function. 

lVj 

Let  W.  = {xc  U.  I q(x)  « V } 1 < i < s,  1 < j < t.  If  W ^ 

1,1  1 J 

then  f I is  a composition  of  two  factorable  functions  and  therefore 

Wi,  j 

factorable.  Since  q(D^)  G ^ ^1,]' 
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Since  locally  (that  is:  on  the  appropriate  sets)  piecewise  factorable 

functions  are  factorable, the  previous  discussion  applies.  If  Z and  T are 
as  in  section  1.4  and  f is  a piecewise  factorable  function  then,  T[f]  is  a 
piecewise  factorable  function.  If  we  carry  out  the  replacement  process  for 
the  factorable  representation  of  each  piece,  we  will  get  factorable 
representation  for  T[  f]  on  each  piece.  We  therefore  arrive  at  the  following 
theorem. 


Theorem  3:  Let  S*  and  T satisfy  assumptions  1-5.  Let  f : D C.  Rn  -♦  R 


be  a piecewise  factorable  function.  Let  U^,  . . . , U C R be  sets 
s i 

is  a factorable  function,  1 < i < s. 


9 D C O U and  f 
~i=l  1 


U. 


Let  f f.  T , 1 < i < s be  factorable  representations  for 

i,r  ’ ' T ’ - - 


i.L 


f . If  we  carry  out  the  replacement  process  as  in  Theorem  2 for  each 
i 

sequence  f.  .,...,  f.  T i<i£s>  then  we  get  a piecewise  factorable 

1,1  1,L. 

A (k-l)-n  Dm • k ~ ( r ,i 

- R • 5 F TT  = Tl f TT  J where 


function  , F : D X R 

U . = U.  X Rn  ‘ (k_1)  . 
i l 


|Ui 


|U,J 


We  call  this  function  T[  f]  . 


In  practice,  as  Examples  1 and  2 show,  5 is  a set  of  analytic 

and  piecewise  analytic  functions.  Also, most  functions  we  use  in  applications 

are  piecewise  factorable  and  therefore  piecewise  analytic.  Of  course, if 

one  wants  to  talk  about  piecewise  analytic  functions,  the  pieces  (that  is 

the  sets  U.)  should  be  "nice"  sets. 
i 

Most  computer  programs  that  arise  in  numerical  applications  compute 
piecewise  factorable  functions.  We  can  always  make  T[  f]  * = f.  Therefore 
if  we  leave  the  decision  statement  and  GO  TO  statement  unaltered  by  the 
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replacement  process,  we  will  not  change  the  path  of  execution.  Since 
locally  this  path  defines  a factorable  function  f,  after  the  replacement 


we  will  get  a factorable  function  which  is  (locally)  T[  f]  . 

1.8.  Iterative  procedure: 

Many  functions  are  computed  iteratively.  The  number  of  iterations 
in  the  computer  program  must  be  finite  of  course.  This  number  however 
can  change  according  to  the  values  of  the  arguments.  The  usual  arrangement 
in  iterative  procedure  is  as  follows:  One  prescribes  a tolerance  e and  sometimes 
an  initial  guess.  The  program  then  proceeds  with  the  iterations  until  the 
change  in  the  function  value,  or  the  estimated  error  is  < e . Once  « is 
fixed,  the  number  of  iterations  as  a function  of  the  arguments  is,  in  most 
cases,  piecewise  constant.  Since  interative  procedure  can  differ  considerably, 
we  cannot  say  what  are  the  precise  conditions  that  make  functions  that 
are  computed  iteratively,  piecewise  factorable.  However  by  careful 
study  of  a particular  problem  at  hand,  in  many  cases,  one  can  show  that 
the  function  actually  computed  is  in  fact  piecewise  factorable. 

In  such  a case  if  one  computes  derivatives  of  that  function, 
one  actually  computes  the  derivatives  of  a piecewise  factorable 
function.  These  derivatives  might  not  be  a good  approximation 
to  the  derivatives  we  had  in  mind.  However  many  times  one  can  use  the 
following  classical  theorem  (see  [ 15]  ). 
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Theorem : 


If  t.  is  a sequence  of  analytic  functions  in  the  complex 

plane  end  f.  — f uniformly  on  a closed  disc  D(Xq,  r)  then  f is  analytic 


{ k ) ( k ) 

in  the  disc,  fv.  -*  f'  , uniformly  and  Vzt  D(x0,r) 

rr  lf(k>(z)  - f'k)(z)  | <-^t  sup  if  (w)  - f(w)  | . 
k-  J rkw,D(x0,r)  J 

1 . 9 Taylor  Series  Expansion  of  Implicit  Functions. 

Let  f « Ck , ( k > 1),  f :RX  R‘n  — Rn  assume  that  f( tQ , xQ)  = 0,  t^  e R, 
x0  t Rn  and  r = f \t  ,Xq)  exists.  Then  the  system  of  equations  f(t,x(t)) 

n 

= 0 defines  implicitly  a unique  function  x:R-»  R , xe  C $ x(tQ)  * x0  an^ 

f(f#x(t))  = 0 in  the  neighborhood  of  ft  .x  ) 

0*  0 

Moreover,  from  Theorem  20.  3 in  [6,  Ch.  5,  p.  329-332]  it  follows  that 


if 


g (t)  = f(t,  l — — (t-t  )J) 

1=0 


then 

i)  for  m < i 


,m 


g.(t) 


dt 


m i 


= 0 


(t0 


i = 1,2,.  . . ,k 


ii)  *U,(t0)  = - r — gi<t) 

dt 


+ There  is  an  error  in  the  statement  of  the  theorem  in  [6]  (Theorem  20.  3, 
Ch.  5,  p.  329). 

.i 


and  not 


xU)  = -r  (g  (t>) 
dt 


(i)  _Lr  JL-  ln  ttw 
x = ‘ i«r  i ))  • 
dt 
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As  a corollary  we  have:  Let 


x^(t 


* .(t-tj 


then 


g.(t)  = f(t,  l -yr*(t-t QV  + v 


,1 

d ( i ) ^ 

•r  7 V = * <V  - * 

dt  '=‘o 

n 


tor  any  vector  x e R . 


The  above  discussion  shows  the  following: 


Theorem  4: 


Let  f be  a factorable  function  (with  X as  in  Example  1). 

0-  V 


f : R X Rn—  Rn.  Assume  that:  ft  CS  (s  > 2),  f(t  x.)  = 0 and 


r = fX(V  V eXiStS> 

Assume  we  know  x(tQ),  [ x(tQ)]  [x(t0)].  (i  < s). 

If  we  take  any  vector  in  Rn  as  an  initial  guess  for  [x(tQ)]i  and  we 

get  the  Taylor  series  expansion  of  the  modified  Newton  method  (r 

considered  as  constant  matrix)  x,  , . = xt  - r f(t  x.  ) then  after 

k+1  k O’  k 

one  iteration  we  will  get  [x  ( )]  exactly: 

Therefore  we  use  the  following  procedure:  We  start  by  the 

regular  Newton  method  to  find  xQ  such  that  f(tQ,  xQ)  = 0.  Then  we  set 
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r 


r f'*(t  x ) and  get  the  Taylor  series  expansion  of  the  modified 
x 0’  0 

Newton  method.  In  the  ith  step  we  compute  derivatives  up  to  order  i, 
taking  0 as  the  initial  guess  for  [x(tf))]..  We  continue  the  same  way 
until  we  get  derivatives  up  to  the  desired  order. 

1. 10  Computation  of  Sparse  Gradients 

a)  Band  gradients:  Many  times  the  gradient  matrix  is  in  band  form  ( for 
example  in  the  numerical  solution  of  a two-point  boundary  value  problem).  In 
such  a case  the  gradient  can  be  computed  with  great  saving  in  time  and 

( 

space. 


Let  f : R"  -*  Rn  oe  a factoraole  function  (piecewise  factorable). 

and  assume  that, for  all  j,f3  depends  only  on  (x1;  where  |i-j|  <k}. 
Assume  that  we  want  to  compute  the  partial  derivatives  of  f with  respect 
to  x£  Rn.  Instead  of  reserving  n+1  words  for  each  variable  and  starting 
with  the  matrix: 


V. 

One  needs  to  reserve  only  2k+2  words  for  each  variable  and  start  with  the  matrix 
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The  implementation  of  automatic  computation  of  partial  derivatives 
of  FORTRAN  functions  (GRADIENT)  described  in  this  report  does  not  use 
the  above  method.  We  plan  to  implement  this  method  in  the  near  future. 

2 . The  Implementation. 

2.1.  Introduction. 

In  earlier  sections  it  was  pointed  out  that  most  computer  functions 
and  subroutines  used  in  numerical  computations  compute  (represent) 
piecewise  factorable  functions.  It  was  shown  that  every  sequence 
representing  a piecewise  factorable  function  can  be  transformed  into 
another  sequence  which  represents  the  original  piecewise  factorable 
function  and  its  total  or  partial  derivatives.  The  translation  process  is 
merely  a replacement  process  and  can  be  carried  out  mechanically. 

In  order  to  implement  such  a replacement  process  one  needs  a 
processor  that  will  do  the  following: 

a)  Break  the  subprogram  into  a sequence  of  one  step  terms. 

b)  Replace  each  term  by  a body  of  code. 

c)  Expand  each  program  variable  into  a vector.  The  size  of  that 
vector  will  depend  on  the  order  of  the  operator  implemented. 

d)  The  processor  should  leave  the  control  structure  unaltered, 
that  is:  Do  loops  and  IF  statements  should  be  left  unaltered. 

There  are  three  principal  ways  to  implement  the  replacement  process. 

a)  Macro  expansion. 

b)  Replacing  each  term  by  a subroutine  call. 

c)  Using  an  interpreter. 


-2t- 


We  chose  to  use  the  second  method  mainly  because  we  could  use  an 


existing  precompiler:  AUGMENT.  We  also  feel  that  the  second  method  is 
the  most  powerful  and  flexible  of  the  three. 

2.2.  The  AUGMENT  Precompiler. 

Since  we  are  using  AUGMENT  as  the  main  tool  for  the  implementation 
of  automatic  differentiation,  we  find  it  appropriate  to  give  a very  short 
description  of  the  function  and  use  of  AUGMENT.  However,  in  order  to 
understand  fully  how  to  use  it,  the  user  should  read  [3]. 

The  AUGMENT  precompiler  was  designed  to  simplify  the  use  of 
nonstandard  data  types  in  FORTRAN.  AUGMENT  enables  one  to  define 
new  data  types  and  operations.  It  enables  one  to  write  FORTRAN  programs 
using  these  new  data  types  as  though  they  were  standard.  AUGMENT 
input  consists  of  programs  written  in  "extended"  FORTRAN,  that  is; 
FORTRAN  programs  using  nonstandard  data  types,  operators,  and  functions. 
AUGMENT  translates  the  input  programs  into  standard  FORTRAN  programs 
with  the  nonstandard  constructs  translated  into  subroutine  and  function 
calls.  The  supporting  package  (that  is,  the  above  subroutines  and 
functions)  implement  the  operations  with  the  nonstandard  data  types. 

In  order  to  implement  a new  data  type  with  AUGMENT  one  has  to  do 
two  things 

1)  write  a package  of  subroutines  to  implement  the  operations  and 
functions  defined  on  the  new  data  type, 

2)  write  a description  deck  which  describes  the  new  data  type. 
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In  the  description  deck  one  provides  AUGMENT  with  the  following 


information: 

a)  The  name(s)  of  the  new  data  type(s). 

b)  The  number  and  type  of  computer  words  to  reserve  for  each 
nonstandard  variable. 

c)  The  set  of  operations  and  "standard"  functions  defined  on  the 
new  data  type. 

d)  The  names  and  calling  sequences  of  the  subroutines  and  functions 
that  implement  the  above  operations. 

e)  The  relations  between  the  new  data  types  and  other  data  types 
(standard  and  nonstandard). 

For  more  detailed  information  see  [3]. 

Z . 3 . GRADIENT  and  TAYLOR  Packages. 

This  report  describes  the  implementation  of  two  types  of  differentiation: 

1)  TAYLOR:  Automatic  Taylor  series  expansion  of  FORTRAN  functions. 

2)  GRADIENT:  Automatic  gradient  computation  of  FORTRAN  functions. 

The  automatic  differentiation  is  implemented  by  providing  two  new  data 
types:  TAYLOR  and  GRADIENT.  The  operations  defined  on  the  new  data 
types  are  the  arithmetic  operations  and  almost  all  the  standard  FORTRAN 
functions.  Each  subroutine  that  implements  a nonstandard  operation 
computes  (represents)  the  corresponding  factorable  function  G of  Theorem  1. 

II 
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Suppose  one  wants  to  compute  the  gradient  of  a function  f,  which 
is  a function  of  n variables.  One  has  only  to  write  a FORTRAN  function 


(subroutine)  that  computes  f.  In  the  function  or  subroutine  code,  one 
declares  the  arguments  and  other  variables  'including  the  function  itself) 
as  type  GRADIENT.  Then  one  submits  this  function  with  the  description 
deck  which  is  part  of  the  package)  to  AUGMENT  as  data.  The  output 
from  AUGMENT  will  be  the  desired  subroutine.  AUGMENT  will  translate 
the  function  into  a FORTRAN  subroutine  written  in  ANSI  standard  FORTRAN, 
declaring  each  GRADIENT  variable  as  a REAL  vector  of  dimension  n + 1 
(and  each  k vector  as  an  (n  + 1)  x k REAL  array  and  so  on).  Each 
arithmetic  operation  or  function  call  will  be  translated  to  a call  to  the 
appropriate  subroutine.  The  translated  subroutine  together  with  the  sub- 
routines provided  by  the  GRADIENT  package  will  compute  the  gradient  of 
the  function  f at  any  desired  point. 

Below  we  give  a detailed  description  of  the  two  packages  and  their 
use.  The  complete  listing  of  the  supporting  packages  and  the  description 
decks  is  given  on  a microfiche  card  at  the  back  of  this  report.  Most  of 
the  details  of  the  two  packages  are  the  same:  Anything  that  is  said  below 
applies  to  both  packages,  unless  the  contrary  is  explicitly  stated.  We 
will  use  the  term  VARIABLE  for  either  type  GRADIENT  or  TAYLOR  and  CONSTANT 
for  types  REAL,  INTEGER  or  DOUBLE  PRECISION.  The  relations  between 
VARIABLE  and  type  COMPLEX  are  undefined. 
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TAYLOR  and  GRADIENT  Variables: 

Each  GRADIENT  variable  is  a REAL  vector  of  dimension  N + 1 where 
N is  the  number  of  the  independent  arguments.  The  first  word  holds  the 
variable  value  and  the  (I  + 1 ) th  word  holds  the  partial  derivative  of 
that  variable  with  respect  to  the  1th  independent  argument. 

Each  TAYLOR  variable  is  a Real  vector  of  dimension  N + 1 where 
N is  the  highest  normalized  derivative  to  be  computed.  The  (I  + l)th 
place  holds  the  Ith  normalized  derivative,  I = 0,  1 , . . . , N. 

Arithmetic  Operations: 

All  arithmetic  operations  between  VARIABLES  and  all  arithmetic 
operations  between  VARIABLE  and  CONSTANT  are  legal  except  INTEGER 
raised  to  a VARIABLE  power.  Since  the  recurrsion  relations  that  replace 
arithmetic  operations  between  CONSTANT  and  VARIABLE  are  simpler  than  the 
general  relations,  separate  routines  are  provided  to  implement  the 
arithmetic  operations  between  CONSTANTS  and  VARIABLES.  Conversion  of 
CONSTANT  to  a vector  format  is  done  if  there  is  a statement  of  the  form 
V - c where  V is  a VARIABLE  and  c is  a REAL  expression,  or  if  there 
is  a reference  to  a conversion  function  (see  Conversion  routines). 

Standard  Functions: 

In  Table  1,  we  list  the  standard  functions  that  are  implemented  in 
the  two  packages.  One  can  easily  add  other  functions  to  that  list  by  adding 
their  names  to  the  description  deck  and  writing  subroutines  to  implement 
them. 
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Table  1 . 


function  reference/  Function  Suffix 


ACOS(x) 


cos  (x) 


Comments 


ALOG(x)' 


fn(x) 


ALOGlO(x) 
AMAXlfx,  y)^ 


AMINl(x,  y)J 


ATAN(x) 


ASIN(x) 


logio(x) 

max(x,y) 


min(x,y) 
tan  \x) 
sin  ^(x) 


function  of  two  arguments 
only 


CBRT(x) 


COS(x) 


COSH(x) 

COTAN(x)1 


cosh(x) 


cotan(x) 


EXP(x) 


exp(x) 


LOG(x)' 


MAX(x,y) 


MIN(x,y) 


max(x,y) 


min(x,y) 


function  of  two  arguments 
only 


SIN(x) 


sin(x) 


SINH(x) 


sinh(  x) 


SQRT(x) 


TAN(x) 


tan(x) 


Not  an  ANSI  standard  function. 


ee  automatic  typing. 


Automatic  Typing: 


This  package  provides  the  same  feature  that  exists  in  many  FORTRAN 
compilers,  automatic  typing  of  functions;  that  is,  the  type  of  function 
used  is  determined  by  its  arguments.  Thus,  for  example,  we  have  defined 
LOG  and  A LOG  to  be  names  of  the  function  which  takes  the  logarithm 
of  an  argument  of  type  VARIABLE. 

Conversion  Functions: 

There  are  three  subroutines  that  implement  conversion  from  CONSTANT 
to  VARIABLE,  one  each  for  types  REAL,  INTEGER  and  DOUBLE  PRECISION. 

These  routines  can  be  referenced  in  the  original  program  by  the  use  of 
the  conversion  function:  CTTYL(-)  in  TAYLOR  and  CTGRD(-)  in  GRADIENT. 

The  function  accepts  all  three  standard  types  as  arguments  (see  automatic 
typing).  Automatic  conversion  is  invoked  only  for  type  REAL.  That  is: 
the  statement  V = const  is  legal  only  if  the  const,  is  of  type  REAL. 

Norm  Function: 

It  is  sometimes  convenient  to  test  the  distance  between  two  VARIABLES 
(for  example  in  a test  of  convergence).  Since  the  relational  operators  compare 
only  the  first  words  of  the  VARIABLES  they  cannot  be  used  for  that  purpose. 

The  packages  provide  a function  NORM  that  computes  the  distance  of  a 
VARIABLE  from  the  0 vector.  In  TAYLOR  package,  the  function  NORM  is  a 
function  of  two  arguments,  TAYLOR  and  REAL 

NORM(v,t)  = max  l[  vj.lt1. 

0 < i < N 1 


r 

In  GRADIENT,  NORM  is  a function  of  one  argument 

NORM(V)  = max  (J V.  |). 

0<i<N  1 

In  both  packages  the  function  is  implemented  as  REAL  function. 

Logical  Statements : 

The  relational  operators  can  be  used  to  compare  two  VARIABLES,  or 
VARIABLE  with  type  REAL.  The  comparison  is  done  between  the  first  words 
of  the  VARIABLES  or  between  the  first  word  of  the  VARIABLE  and  type  REAL. 

The  comparison  operators  are  implemented  as  LOGICAL  functions. 

Other  Subroutines: 

The  packages  provide  two  additional  subroutines: 

1)  Error  handling  subroutine  (see  our  later  discussion  of  Error  handling). 

2)  Copy  subroutine. 

The  copy  subroutine  implements  the  statement  A = B,  A and  B VARIABLES. 
Subroutine  Names: 

The  names  of  all  subroutines  in  both  supporting  packages  are 
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The  suffix  of  a subroutine's  name  depends  on  the  function  or 
operation  the  supporting  routine  implements.  In  Table  1,  we  give  the 
suffixes  of  the  routines  implementing  the  "Standard"  functions.  The 


suffix  of  a routine  implementing  arithmetic  operations  is  given  systematically. 
The  first  letter  describes  the  operator.  A for  +,  S for  M for  *f 
D for  /,  and  E for  **.  The  next  two  letters  describe  the  operands: 
first  the  left  operand  and  then  the  right  one.  The  letter  R stands  for 
REAL,  I for  INTEGER,  D for  DOUBLE  PRECISION  and  V for  VARIABLE.  So 
MW  will  be  the  suffix  of  a routine  that  implements  (VARIABLE)  * (VARIABLE) 
and  DDV  of  a routine  that  implements  (DOUBLE  PRECISION)/( VARIABLE) . 

The  suffix  of  the  LOGICAL  functions  implementing  the  relational 
operators  is  composed  of  the  two  letters  representing  the  operator  and  the 
letter  V.  So  . LT.  is  implemented  by  a function  with  suffix  LTV.  The 
suffix  of  the  routine  implementing  the  norm  function  is  NRM,  Error  routine  - ERR 
and  copy  routine  - CPY. 

2. -4.  Using  the  Package  with  AUGMENT. 

Writing  the  Source  Code: 

In  order  to  get  derivatives  of  a function,  say  the  Taylor  series 
expansion  of  a function  f,  the  user  should  write  an  "extended"  FORTRAN 
function  or  subroutine  that  computes  f.  All  legal  FORTRAN  constructs  can 
be  used.  All  program  variables  which  depend  on  the  independent  variable 
should  be  declared  as  type  TAYLOR,  including  the  function  itself.  Program 
variables  which  do  not  depend  on  the  independent  variable  can  be  of  any 
other  type. 
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External  functions  and  subroutines  can  be  used  as  part  of  the 


computation.  Functions  must  be  declared  as  type  TAYLOR.  External 
functions  and  subroutines  can  be  translated  separately.  However,  one 
has  to  make  sure  that  the  number  of  computer  words  reserved  for  each 
TAYLOR  variable  in  the  external  subroutine  (function)  is  the  same  as 
the  number  of  words  reserved  in  the  calling  routine. 

The  Description  Deck: 

The  next  step  is  to  submit  the  source  deck  with  the  description  deck 
as  data  to  AUGMENT.  See  Appendix  C for  the  deck  structure.  The 
description  deck  is  supplied  with  the  package. 

However,  since  the  number  of  words  reserved  for  each  VARIABLE 
changes  from  problem  to  problem,  this  number  has  to  be  put  into  the 
description  deck  each  time.  To  make  it  easier,  the  description  deck  was 
split  into  two  parts:  HEAD  and  BODY.  The  number  of  words  to  reserve  is 
inserted  in  between  the  two  parts.  This  number  should  be  an  integer 
constant.  Column  1 in  the  card  holding  this  number  must  be  left  blank. 

Order  of  Differentiation: 

The  routines  in  both  packages  were  designed  to  implement  any  "order" 
of  differentiation  without  the  need  to  be  recompiled  every  time  the  "order" 
is  changed.  The  "order"  of  differentiation  is  provided  to  the  package  through 
a common  block.  In  TAYLOR  by  COMMON/DEGREE/N  and  in  GRADIENT  by 
COMMON/ORDER/N  N is  the  order  of  the  operator,  that  is:  if  N = 1 
the  package  will  compute  function  values  only;  if  N = 2,  function  values 
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and  first  derivatives  (or  first  partial  with  respect  to  one  variable);  and 
so  on.  The  routines  in  the  package  do  not  check  that  there  is  enough 
space  provided  for  the  VARIABLES.  However  there  is  a check  that  N > 1. 

The  order  can  be  changed  at  run  time  but  care  should  be  taken  not  to 
exceed  the  number  of  words  provided  for  each  VARIABLE. 

Working  Space: 

Some  routines  in  TAYLOR  need  work  space  and,  since  they  are 
designed  to  handle  any  order  of  differentiation,  the  work  space  has  to  be 
provided  by  the  user.  The  work  space  is  provided  through  four  common  blocks. 
COMMON/WORKl/WORKl(N) 

COM  M0N/WORK2/WO  RK2  ( N) 

COMM0N/WORK3/WORK3(N) 

CO  M MON/WORK4/WO  RK4(  N) 

N should  be  the  highest  order  of  differentiation  used.  The  GRADIENT 
package  does  not  require  work  space. 

Using  the  Translated  Routine: 

The  translated  routine  is  a FORTRAN  subroutine  that  gets  as  input 

the  value  of  its  arguments  and  their  derivatives,  and  gives  as  output  the 

value  of  the  function  and  its  derivatives.  For  example,  in  the  Taylor 

package,  if  t is  the  independent  variable,  then  t,  1,0,0...  is  the 

Taylor  series  expansion  of  t.  In  the  Gradient  package,  if  x is  the  j* 

independent  variable  then  =1  and  = 0 for  i 4-  j. 

ax.  dx, 

J i 
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In  this  example  we  compute  the  first  9 terms  of  the  Taylor  series 


expansion  of  f(t) 


exp(cos(-4t))  + aictanfsin 


l 


(t)) 


at  the  point 


. 5. 


a)  The  source  code: 


TAri  or-.  ruNcrwM  funct) 

TAYLOR  T 

FUN-  EXE’  ( COS  ( 4 . #T  ) ) + AT  AN  (SIN  ( T ) **2  ) 

RETURN 

F.  ND 


b)  The  translated  code: 

SUBROUTINE  TUN  (T,  TYLRlET ) 


c 

PROC!  Sum  BY  AUGMENT t VERSION 

c 

-• 

TEMPORARY  STORAGE  LOCATIONS  ~ 

c 

TAYLOR 

REAL 

TYLTMP<9»2) 

c 

-■ 

LOCA!  VARIABLES  

c: 

TAYLOR 

REAL 

TYERCSC9) 

c 

. 

GLOBAL  VARIABLES  

c 

TAYl OR 

RE  Al 

T ( 9 ) ? 

I Y L Rl  T ( 9 > 

c 

■ -----  TRANSL  ATED  PROGRAM 

CALL 

TYL  MRU 

(4.  ,T?TYLTMP(I  , .1  ) ) 

CAT  r 

TYLCOS 

( TYL IMP ( I v 1 > r T YLTMP < 1 r 1 > ) 

CA1  1 

TYL  TXE 

( TYLTMP ( 1 » 1 ) t TYLTMP ( 1 , 1 ) ) 

CALI 

TYL SIN 

( T » TYL TUP (1,2)) 

CA!  I 

TYLEUT 

< TYLTMP  (1,2)  ,2,  TYL  TMP  (02)) 

CAl  1 

TYLATN 

< TYL  TMP  C L , 2 > * TYI  TMP (1,2) ) 

CAl  1 

TYI  AVV 

( TYLTMP < 1 . L > . TYLTMP (1 ,2) , TYL RES) 

GO  TO  30000 

r 

... 

RETURN  CODE  

30000 

CONTINUE 

CAl  1 

TYLCF'Y 

< TYLRES , TYLRLT) 

El  TURN 
END 


« 


13  li'"' 11  - 

S?.;J 
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c)  Main  program: 


I' T MIMS  ON  T(9).nif’'  9) 

( (UUU'rl  /HI  Pfvf  ! / N 
COMMON  /UORM/W1  (9) 

•COMMON  /UflRI  ?/W?  ( v ) 

COMMON  /U0RN3/U3' V > 
rOMNON  /UORI 4/W4  ( 9 ) 

N 9 

T ( ? ) ■ 1 . 

00  10  I 3,9 
10  1(1)  ;0  . 

COL!  rUN(1,rUNC) 

PR1NT91  , T(  1 > , (FUNCML  ) ,L=1 ,9) 

91  FORMAT (IX#  'CXAMPIFJ  TAYLOR  SERIES  EXPANSION' ,/,3X» 'T=' >FA.3//» 

* (3X»F13.5>> 

STOP 

END 


d)  Output: 


EXAMPLE!  TATI  OR  SERIES  EXPANSION 
T=  .000 

.osr.t. i ioo 

-.  1599!)  i 01 
. A9r:;i  i oi 
7943 |01 
- .342VA-* 01 
.?r,40A10? 

- . 59899+0? 

. ?8534 ! O? 

. 1076?+03 


In  Appendix  B we  give  a more  complicated  example. 

Error  Handling: 

In  general  the  packages  do  not  check  that  the  arguments  are  within 
the  domain  of  definition  of  the  functions.  Checking  is  done  only  for  division 
by  REAL,  INTEGER  or  VARIABLE  types.  The  packages  also  provide  the 
capability  to  specify  what  constitutes  division  by  zero.  If  the  absolute  value 
of  an  argument  to  division  routine  is  smaller  than  or  equal  to  specified 
value,  error  occurs.  This  value  is  set  initially  to  0.0  by  a DATA  state- 


ment. However  it  can  be  changed  in  runtime  through  common  block 


-ZLRO/EPS.  The  routines  in  the  package  also  check  that  the  order  of 
differentiation  is  at  least  1.  In  case  an  error  is  discovered,  the  error 


routine  is  called  (suffix  ERR).  The  error  routine  prints  error  messages 
and  periorms  a walk  back"  (which  traces  the  sequence  of  subprogram 
calls  back  to  the  main  program)  and  stops. 

Mon  ANSI  FORTRAN  Parts: 

No  special  care  was  taken  to  comply  with  some  of  the  restrictions 
imposed  by  ANSI  FORTRAN,  for  two  reasons: 

a)  Some  of  the  features  in  the  packages  could  not  have  been 
implemented  otherwise. 

b)  Some  of  the  restrictions  violated  seem  to  us  arbitrary  and  unreason- 
able . 

Moreover  they  do  not  exist  in  most  production  compilers. 

Below  we  give  the  list  of  non  ANSI  FORTRAN  constructions  used. 

1)  The  set  of  ’ standard"  functions  used  is  larger  than  the  set  in 

ANSI  FORTRAN.  Functions  which  are  not  in  the  Standard  are  flagged 
in  the  function  table  by  f.  The  corresponding  routines  could  be 
deleted  or  modified  to  fit  different  systems. 

2)  Tne  work  space  to  TAYLOR  is  provided  through  common  blocks 
and  therefore  the  common  block  sizes  in  the  TAYLOR  routines 
would  be  different  from  the  block  sizes  in  the  main  program. 

3)  The  walk  back  routine  used  in  the  error  handling  routine  is  non- 
standard and  has  to  be  changed  in  other  systems. 
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4)  The  REAL  variable  EPS  appears  in  common  block  /ZERO/  and  in 


expressions  as  subscripts  to  arrays. 

o)  The  function  ATAN2  is  not  implemented. 

7)  No  care  was  taken  not  to  mix  INTEGER  with  REAL. 

Using  the  Package;  Summary: 

Assume  one  has  all  the  package  subroutines  in  relocatable  form,  so 
they  can  be  used  as  library  routines.  In  order  to  get  the  derivatives  of  a 
function  f one  has  to  do  the  following: 

A.)  Write  a FORTRAN  subroutine'  (function)  that  computes  the  function 
f.  In  the  subroutine,  declare  all  FORTRAN  variables  that  depend  on  the 
independent  variables  as  type  TAYLOR  (or  GRADIENT).  The  rest  of  the 
variables  can  be  of  any  other  type. 

B)  Insert  the  number  of  computer  words  reserved  for  each  VARIABLE 
into  the  description  deck. 

C;  Submit  the  source  deck  with  the  appropriate  description  deck  as 
data  to  AUGMENT  (See  Appendix  C and  also  [3]). 

D)  Submit  the  output  from  AUGMENT  as  data  to  the  FORTRAN  compiler. 

E)  Call  the  subroutine  with  the  desired  arguments. 


+ One  can  use  main  programs  too  but  that  is  a more  complicated  way  of  doing  it. 
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Remarks 

l)  Always  remember  that  the  initial  values  of  the  derivatives  of  the 

input  variables  have  to  be  provided  too. 

1)  All  the  restrictions  that  apply  to  the  use  of  AUGMENT  apply  to  the 

packages . 

3)  AUGMENT  does  not  translate  I/O  and  DATA  statements.  Their 
translation  has  to  be  done  by  hand. 

4)  This  report  is  by  no  means  a substitute  for  AUGMENT  user  informa- 
tion manual  ( MRC  TSR  #1469  [ 3]).  The  user  should  be  familiar  with  that 

report. 
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Appendix  A 


Recurrence  Relations  for  Taylor  Series  Expansion: 


a) 


Z = X ± Y 


b) 


c) 


d) 


f) 


[Z]j  = [X]j  ± [Y]J 


Z = X • Y 

[zij  = h 

J k =0 
Z = X/Y 


[XJk  • [Y] 


I-k 


J = 0,1,... 


J = 0,1,.. . 


J-l 


fZl  = l/Y  { [X]  - £ [Z]  • [Y]  } J = 0,1, 

1 1 k = 0 k 


Z = X 


1)  I > 0 , I = V Lq  2' 

s=,°  s 

n L • 2 

[z]T  - [ I I xs  ] 

J s=o  J 

2)  I =0  ; Z = 1 


I integer 
S 


Lg«  {0,1} 


J = 0,1, 


3)  I < 0 [Z]  = [1/X_I]J 


e)  Z = X 


J = 0,1,... 

a real  constant 


1-1 


[Zb  = 1/X  • £ ((a( j-k)-k)/j)  • [Z],  • [X] 

J k=0  J 


Z = X 


Y 


J = 1,2, 


[Z)]  = fEXP( Y • LOG(X))l  j = 0,1,.  . . 


g)  Z = LOG(X) 


[Z\  = [X]j/X 


1-1 


[Z]  = 1/X  {[X]  - £ (( J-k)/j)[X]  • [Z]  J = 2,3,.. 

J J k=0  J k 

h)  Z = EXP(X) 

J-l 

[Z]  = V ( ( J- k )/T ) [Z],  • [X]  J = 1,2,... 

J k =0  J'K 
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1) 


j) 


k) 


n 


Z = SIN(X),  W = COS(X) 

J-l 

[Z]  = V (( j-k)/j)[w]  .[XI 

J k=0  J'K 

J-l 

[W]  = - V (( J-k)/j)(Z]  • [X] 

J k"o  k J-k 

Z = SINH(X),  W = COSH(X) 

J-l 

[ZL  = J <(J-k)/j)[W]  ■ [XI 

1 k^O  k j'k 

J-l 

[W]  = V (( J-k)/j)  • [Z]  • [X] 

J k=0 

Z = ATAN(X) 

Let  V = X2  + 1 W = l/V 

[Z]j  = 2/(J-k)/j)  [W]k  • [X]J  k 
Z = ASIN(X),  Y = ACOS(X) 

Let  V = 1 - X2  W = 1 /\TT 

J-l 

[Z ] = X (( J-k)/j)[Wl  • [X] 

1 k=0  k J 

J-l 

[Y]t  = -£  (( J-k)/j)  [W]  • [X]  , 

J k"o  k J'k 

TAN(X)  = SIN(X)/ COS(X) 

\[T  = X a 

TANH(X)  = SINH(X)/ COSH(X) 


J = 1,2, 


J = 1,2,. . . 


J = 1,2,.. . 
J = 1,2,... 
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Appendix  B 


Example 

The  following  example  was  taken  from  a homework  assignment  given  in 
an  optimization  class  at  the  University  of  Wisconsin.  This  example  is  simple 
but  quite  typical  of  problems  that  arise  in  applications.  We  have  to  compute 
the  gradient  of  a function  which  can  be  easily  described  by  a computer  program, 
but  whose  explicit  expression  is  quite  hard  to  obtain.  When  one  tries  to  com- 
pute the  gradient  numerically,  one  runs  into  convergence  problems.  In  this  ex- 
ample we  show  how  easy  it  is  to  get  the  gradient  of  such  a function  by  using  the 
GRADIENT  package. 

The  Problem 

We  have  a missile  on  the  north  pole  of  a ball  of  radius  1 and  we  want 
to  fly  to  the  south  pole.  (The  units  are  chosen  in  such  a way  that  all  the  con- 
stants are  1).  Because  the  problem  is  symmetric  we  only  have  to  solve  a two 
dimensional  problem. 


The  motion  equations  are: 


A 

dt2 


+ u 


r = (x,y) 


Let  t = ar  then 


d 

dt 


(*\ 

I y 


= a 


t 

n 

-X 


, 2 2.3/2 

\(x  + y ) y 


af(x,y,|,r1,u1,u2) 
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Discretize  by  the  Euler  method 


Finally,  solve: 

min  h(a,u^0,u2}0,u1>8,u2>8,ulj9,u2>9)  . 

= 12°[<X10)2  + (y!0  + 1)2  + (410)2  + (T110,2t  + (U12,  0)2  + (U2,0)2 

9 

20 


+ ,Ul,8|2  + ,U2,8|2  + |Ul,9|2  + ,Ui,9,2+  + 2.2 

X5  Y 5 


} • 


All  u.  ,s0  u . = 0 for  1 < j < 8 . 

1,1  2,J 

Use  the  variable  metric  algorithm  to  solve  the  above  minimization  problem. 
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1. 


INPUT  DECK 


SX3T  MRC*LI3. AUGMENT 
SA33  GK«OIFF.OESC-GRD/HEAO 
3 

iA  33  GK.DIFF.DCSC-GRD/SODY 
•3EGIN 

IFOR. IS  TEST1 

IMPLICIT  GRADIENT  (A-H.0-2) 

GRADIENT  FUNCTION  FL TFN ( ALF A . UO .U3. U9 ) 

DIMENSION  U0(2).  U 3(  2 1 . U3  < 2 I .U I 2 . 1C  I . XI 11 1 • Y 1 11 1 • VY 1 11 ) t VX 1 11 1 

X(ll=C.O 

Ylil=1.0 

V X( 1 1 = C.Q 

V Y( 1 1 =C.Q 
Hz. I 

DC13I=1.1C 
U(l. 11=0.0 
U ( 2 • I 1:3.3 
1C  CONTINUE 

Ulltl IzUOIll 
U (2 .1 1 =UC  (2  ) 

U(l, 31=03(1) 

U(2 .3 IzUS 12  ) 

U(I. 10)=U3( I) 

UI2.10IZU3I2I 
3020  1=1. 1C 

RSzlX (I>»  »2*Y II I «»2I «»1.5 
XII*1 1=X( II *H»ALFA*VXII I 
YtI*il=Y(Il*H»ALFA»VY(H 
VX(I*Ilrv XIII *H*ALFA»IUI1.I)-Xtll /RSI 
V Y( I*11=VY( II *H*ALFA»IUI2*I1-YI II/RS1 
20  CONTINUE 

FLTFN  = 20. •<  XI 11 1 »»2* ( Y( 11 l*li»»2*VX 111) • • 2*  VY  1 1 1 1 • • 2 1 

* ♦ UC(1 ) •• 2*U0 (2) ..2+U8 1 1) ••2*U8 <2 !• »2*U9( II **2*U3 121 »»2 

* ♦ (ALFA.  .2  1/13.0  *20 .0  / 1 XI  6 1 * * 2*  Yl  6 1 **2 1 
RETURN 

ENO 

• END 
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Remarks. 


1)  @XQT  MRC*LIB.  AUGMENT  starts  the  execution  of  AUGMENT. 

2)  (gftDD  GK*DIFF.  DESC-GRD/HEAD,  adds  the  card  image  of  the  HEAD  part  of 
the  description  deck  into  the  run  stream. 

Similarly  (§ADD  GK*DIFF.  DESC-GRD/BODY  Adds  the  BODY  part. 

3)  The  function  is  translated  into  a subroutine.  The  function  and  gradient  values 
are  stored  in  the  last  argument  of  the  subroutine. 

4)  The  translated  subroutine  can  be  used  to  compute  function  and  gradient  values 
or  function  values  alone.  (See  Order  of  Differentiation) 

5)  The  rest  of  the  computation  details  are  omitted. 
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OUTPUT  FROM  AUGMENT 


SUBROUTINE  FLTFN  t ALT  A . UC  .U  3.  U9  • GRORLT) 

C :::::  PROCESSED  ST  AUGMENT.  VERSION  AH  ::::= 

C TEMPORARY  STORAGE  LOCATIONS  

C GRADIENT 

REAL  GRDTMPI3.3) 

C LOCAL  VARIABLES  

INTEGER  I 

C GRADIENT 

REAL  H ( 3 ) * RS  ( 8 ) . U13.2.1C).  VXt8.ll).  VYtS.lU.  X t 3 , 11 ) . YI8.11)* 

• GRDRCSI 3) 

C GLOBAL  VARIABLES  

C GRADIENT 

REAL  ALF  A < 6 ) . U0I3.2).  U316.2),  U9I8.2).  GRDRLT { 8 ) 

C TRANSLATED  PROGRAM 


CALL  CROC  FR  (D.C.X(l.l)  ) 

CALL  GRDCFR  tl.O.Ytl.l)) 

CALL  GRDCFR  {C.C.VX11.U) 

CALL  GRDCFR  I 0. 0. V Y I 1 . 1 > ) 

CALL  GRDCFR  t .1 .HI 

DO  ID  1:1.10 

CALL  GRDCFR  ( 0. G . U t 1 . 1 . I ) ) 

CALL  GRDCFR  ( 0. D . U (1 . 2. I ) ) 

10  CONTINUE 

call  grdcpy  (uca.i)  .ud.i.m 

CALL  GRDCPY  l UG  ( 1 . 2)  . Ut  1 . 2 . 1 ) ) 

CALL  GRDCPY  ( U3 t 1 . 1) . U( 1 . 1 , 9 )) 

CALL  GRDCPY  t U3 11.2) .Ut 1.2, 9) > 

CALL  GRDCPY  C U3 t 1. 1) . Ut 1.1. 1C) ) 

CALL  GRDCPY  t U9 t 1.2) . Ut 1.2, 10) ) 

DO  20  1:1, ID 

CALL  GRDEVI  I X t 1 . 1 ) . 2 . G RD TM P f 1 , 1 ) ) 

CALL  GRDEVI  t Y t 1 . I ) . 2 . GRD TM P ( 1 , 2 ) ) 

CALL  GRDAVV  t GR DT MPt 1 , 1 ) , GRDTMP t 1 .2 ) . GRDTMP t 1 .2 ) ) 
CALL  GRDEVR  t GR  DT  MPt  1 . 2 ' - 1 . 5 . RS  ) 

CALL  GRDMVV  l H.  ALF  A.  G RD  i .IP  l 1 , 1 ) ) 

CALL  GRDMVV  t GRDTMPt 1 . 1 ) . VX t 1 . I ) . GR DTMP t 1 .1 ) ) 

CALL  GRDAVV  t X t 1 , I ) . GRD TMP ( 1 , 1 ) , X 1 1 .1  + 1 ) ) 

CALL  GRDMVV  ( H. ALF A. GRD TM P t 1 , 1 ) ) 

CALL  GRDMVV  t GR DTMP f 1 . 1 ) . VY l 1 , I ) . GR DTMP t 1 , 1 ) ) 

CALL  GRDAVV  t Y ( 1 . I ) . GRD TM P ( 1 , 1 ) , Y 1 1 .1*1 ) ) 

CALL  GRDMVV  l H. ALF A, GRD TMP l 1 , 1 ) ) 

CALL  GRDDVV  IXtl.I). PS, GRDTMP (1.2)  ) 

CALL  GRDSVV  t U 1 1. 1 • I > » G RD  TM  P 1 1 . 2 ) .CRD  TMP  1 1. 2 I ) 
CALL  GRDMVV  l GR DT MP t 1 . 1 ) , GR DTMP ( 1 .2 ) . GRDT MP ( 1 .2 ) > 
CALL  GRDAVV  t VX  1 1 . I)  . GR  DT  MP  1 1 . 2 ) . VX  1 1 .1*1  ) ) 

CALL  GRDMVV  t H. ALF A. GRD TM P l 1 , 1 1 ) 

CALL  GRDDVV  t Yt 1. I ) . RS. GRDT MP t 1 .2 ) ) 

CALL  GRDSVV  t J( 1.2.1) .CRDTMPt 1.2) *GRDTMPtl,2) ) 
CALL  GRDMVV  I GR DT MPt 1 , 1 ) , GR DT MP 1 1 ,2 ) , GROT MP 1 1 .2  ) ) 
CALL  GRDAVV  t VY t 1 . I ) . GR DT MP l 1 . 2 ) . VY t 1 . 1*1 ) ) 
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COPY  AVAILABLE  TO  DOC  WES  NOT 
PERMIT  FULLY  LEGIBLE  FSSSUCTIffli 


continue 

CALL 

GR2E VI 

CALL 

CR2AVI 

CALL 

GR3C v: 

CALL 

GRDA  VV 

CALL 

GRDE VI 

CALL 

GRDA VV 

CALL 

GROG  VI 

CALL 

CR2AVV 

CALL 

GROMRV 

CALL 

CR3C  VI 

CALL 

GRDA  VV 

CALL 

CRGE  VI 

CALL 

GRDA  VV 

CALL 

CRDEVI 

CALL 

GRDAVV 

CALL 

CROC  VI 

CALL 

CRD  AW 

CALL 

GRDCVI 

CALL 

GRDAVV 

CALL 

GRDE VI 

CALL 

GRDAVV 

CALL 

GRDCVI 

CALL 

GRDDVR 

CALL 

GRDAVV 

CALL 

GRDC VI 

CALL 

GRDCVI 

CALL 

GRDAVV 

call 

GRDDRV 

CALL 

CRDAVV 

CO  TO  3000C 

c 

3CGQ0  CONTINUE 

CALL  CRDCPY 

RETURN 

ENO 


1 XI 1,11) ,2,GRDTMP(1 ,1) J 
(Y(l.ll)»l.CRDTMPd,2>) 

(GRDTMP(I.2) . 2,  GROT  HP ( 1 .2  > ) 
t GR3TMP<  1 ,1 ) , GRDTMP  (1  ,2  > . GRDTMP  1 1 .2  > ) 

< VX  <1.11 ) .2 tGRG  TMP( It 1)  I 

( GRDTM?(  1.2 ) .GRDTMP (1 .1  ), GRDTMP (1,1  ) ) 

< VY  d.  11  ) » 2 • GR3 TMPd. 21  > 

( GRDTMPt 1,1) .GRDTMP (1.2 1 .GRDTMP (1.2  I ) 
(20. . G RCTMP (1.2). SPOT MP (1.21) 

( UO 1 1. 1) .2.3RDTMP(1 .1)1 
( GRDTMP  ( 1.2  ) .GRDTMP  (1 ,1  ) . GRDTMPd.il  1 
( U0d.2l  .2.  GROTMPd  .2)  ) 

(GRDTMP(l.l), GRDTMP (1.2). GRDTMP (1.2)1 
{ US (1. 1) .2. GRDTMPI1 .1) ) 

I GROT HT ( 1.2  I .GRDTMP (1 .1 ) .GROT HP (1.1)  ) 
(U811.2) .2. GRDTMP (1.2) ) 

(GRDTHF l 1.1). GRDTMP (1.2). GRDTMP (1,2) ) 
( U3  (1.1)  .2.  GRDTMP  d.l)  ) 

( GR CTMril. 2 ). GRDTMP (1.1), GRDTMP (1.1)  ) 
(US(1.2) .2, GROT M° (1.2) ) 

( GROT MP(  1.1 1 .GRDTMP (1,2 ). GROT  HP (1.2 ) ) 
( ALFA. 2. CRD  TMP( 1,11) 

( GRDTMPI 1 . 1 ) . 10. 0.  GRDTMP  d.l) ) 

( GRDTMP 1 1.2 ) .GRDTMP (1,1), GRDTMP  (1.1) ) 
( X{ 1,5 ).2.CRDTMP(1, 2) ) 

( Y ( 1.6  1.2. GROTMPd.  3) ) 

I GRDTMP< 1.2) .GRDTMP (1.3). GRDTMP (1.3) ) 
( 20.C.SRDTMP(1> 3). GRDTMP (1.3) ) 
(CRDTMPd.l ) .GRDTMP  (1.3).  GRDRES) 

RETURN  CODE  

(GRDRES. GRDRLT) 


COPY  aymuile  TO  OQC  does  not 

PEBWT  FULLY  l£6IBl£  PRODUCTION 
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Appendix  C 


Deck  Arrangement 

The  data  deck  read  by  AUGMENT  has  the  structure  shoun  in  the 
following  diagram: 


/•END  I 


/ Source  Deck  /I 

/ / I 

/ I I 

I * 

I 

/•BEGIN-  I 

I I 

I 

/ Description  Deck  / I 

/ / I 

/ I I 

I I 

I I 


At  the  conclusion  of  processing,  the  translated  program  decks 
are  in  the  output  file  in  80  column  card  image  format. 
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